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By means of molecular dynamics simulations we examine the aging process of a strong glass 
former, a silica melt modeled by the BKS potential. The system is quenched from a temperature 
above to one below the critical temperature, and the potential energy and the scattering function 
C(t w , t + t w ) for various waiting times t w after the quench are measured. We find that both 
qualitatively and quantitatively the results agree well with the ones found in similar simulations 
of a fragile glass former, a Lennard- Jones liquid. 



§1. Introduction 

In the recent years aging effects have become a very 
active field of research in the study of non-equilibrium 
systems. Common features have been found in the aging 
behavior of spin and structural glasses in experiments 
and simulations, and similarities in the theories worked 
out indicate the possibility of a connection between the 
slow dynamics in both types of glasses .El 1 

Even though silica (SiC^), a strong glass former, has 
been examined frequently both in theories and exper- 
iments, comparatively few simulations have been per- 
formed so far on the dynamics of systems like this that 
form an open network structure. In this paper we present 
several results of aging simulations for silica and compare 
them to those obtained for a model for a fragile glass for- 
mer, a binary Lennard-Jones mixture. 

§2. Simulation 

2.1 BKS potential 

The model we used in our simulations is given by the 
potential proposed by van Beest, Kramer and van Santen 
(BKS)© 



U 3 k(r 3k ) = 



Aj k exp(-B jk rj k ) ~ ^jp 
r jk 



with 7 = e 2 /(4-7T£o); the values of the parameters, taken 
from the original publication, are listed in tables I and 
II. As has been shown in previous calculations and 
simulations, this potential reproduces many of both the 
structuralHf and dynamical properties of silica and has 
frequently been used for simulations of the amorphous 
phase. Moreover, it has the advantage of consisting 
of two-body interactions only, so that it can be imple- 
mented in a simulation more efficiently than some other 
potentials proposed for silica that also contain three- 
body terms. Considering that silica forms an open net- 
work where the silicon atoms are the centers of tetra- 
hedra and the oxygen atoms work as bridges between 
them, it is not obvious at all that a potential without 
three-body interactions is suitable for a good description 
of the properties; it has been found, however, that the 



competing two-body forces indeed mimic the three-body 
ones, causing the formation of a network. i3 

Previous simulations^ have shown that the critical 
temperature for this model is T c = 3 330 K, which is a 
good approximation forJhe value found in experiments 
with real silica, 3 221 K.B 



Table I. Interaction constants in the BKS potential. 
j-k A jk [eV] B jk [A" 1 ] C jk [eV . 



O-O 
O-Si 
Si-Si 



o 

Si 



1388.7730 
18003.7572 




2.76000 
4.87318 




175.0000 
133.5381 




Table II. Charges and masses of the ions. 

1j m j M 



-1.2 
+2.4 



15.9940 
28.0855 



For the 0-0 and O-Si interactions, the BKS potential 
has the drawback that it shows an unphysical behavior 
for small distances by diverging towards minus infinity. 
To correct this, it was replaced by a harmonic potential 
for distances below the location of the potential maxi- 
mum (1.4387 A for O-O, 1.1936 A for O-Si). We_con- 
firmed the observations made in earlier simulationsliiP al- 
ready: because of the potential barrier the particles have 
to overcome in order to reach them, these small distances 
only occur in very few particle pairs, so the substitution 
of the potential is not likely to have a relevant effect on 
the quantitative results of the simulation. 

2.2 Parameter choices 

We performed molecular dynamics simulations in a 
cubic system of length L = 25.11 A containing 1089 
(363+726) particles, so that the density was fixed at 

2.3 g/cm 3 . To integrate the equations of motion, the 
velocity form of the Verlet algorithm with time steps of 
1.6 fs was used. The non-Coulombic part of the potential 
was truncated and shifted at a distance of 5.5 A. The 
Coulomb term was calculated by means of the Ewald 
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summation^ splitting it into a sum in real space, 



j<k 



(where erfc denotes the complementary error function, 
erfc(x) = J°° exp(— t 2 ) dt), and one in Fourier 

space, 



2ttL *-f n 2 V 



a 2 L 2 



Qj ex P -T nr : 



(n = |n|). We chose a = 6.5/L; the real-space sum was 
truncated and shifted at 9 A, the Fourier-space one at 
\n\ = 5. Further, details about the simulation can be 
found elsewhere £3 

2.3 Preparation and measurements 

The simulations were done in the same way as previ- 
ous ones for a Lennard-Jones potential:otf a system is 
first equilibrated at a temperature above the critical one, 
then at t — suddenly "quenched" to one below; there it 
is allowed to relax for various waiting times before mea- 
surements of the potential energy and the generalized 
scattering function 

C(t w ,t w +t) = i ^exp(iq • [rj(t + t w ) - 7j(t w )]). 

3 

(where N is the number of particles and q is a wave 
vector) are started, so that examples for both an equi- 
librium constant and a two-time correlation function are 
examined. 

To improve the statistics of the results, we averaged 
over three independent initial configurations consisting 
of randomly distributed particles. Each of those was used 
to produce equilibrium configurations at 6 000, 7 000 and 
8 000 K by first rescaling the particle velocities periodi- 
cally according to the desired temperature and then let- 
ting the system evolve freely for 55 000 time steps. 

Each of the resulting configurations was then used 
for measurement runs at 2 700, 3 000 and 3 200 K: at 
t = 0, the system was brought to the final tempera- 
ture by rescaling the velocities, and the temperature was 
controlled throughout the run by repeating this every 50 
time steps. After waiting 0, 10, 100, 500, 1000, 3 000, 
10000 and 30 000 steps, respectively, measurements of 
the potential energy and the generalized scattering func- 
tion were started and performed every ten steps; the 
scattering function was calculated by averaging over fifty 
randomly chosen wave vectors of absolute value 1.7 A -1 , 
theJirst sharp maximum of the structure factor for sil- 
icarP These runs were executed on eight processors of 
a Cray T3E with a time limit of four hours per run, 
which allowed between 60 000 and 75 000 time steps; 
with the chosen step size, this corresponds to 96-120 
ps (10 _12 sec.) in reality. 

§3. Results 

3.1 Potential energy 

The potential energy shows a fast decay only for a 
short time after the quench; after that it appears to re- 



main constant at a first glance in a linear plot (an ex- 
ample is shown in figure la), but a closer examination 
by means of a logarithmic plot (figure lb) reveals that it 
still decays slowly. 



Fig. 1. Potential energy of the system (shifted by 13 900 eV) as 
a function of the time after the quench from 7 000 to 3 000 K, 
plotted (a) linearly and (b, dotted) logarithmically. The dashed 
line in (b) shows the approximation for long times using a power- 
law dependence with an exponent of -0.22. 
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It is possible to approximate the decay for long times 
by a power law. The exponents for this approximation 
are listed in table III; the average value for all quenches 
is -0.20. 



Table III. Estimated exponents for the time dependence of the 
potential energy for the various combinations of initial (rows) 
and final (columns) temperatures. 





2 700 K 


3 000 K 


3 200 K 


6 000 K 


-0.22 


-0.18 


-0.15 


7000 K 


-0.26 


-0.22 


-0.18 


8 000 K 


-0.28 


-0.18 


-0.16 



Next we compare these results to those obtained in 
analogous simulations for a binary Lennard-Jones (LJ) 
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system. This model glass consists of a binary mixture of 
spheres which interact via an LJ potential with a strong 
r 12 repulsive and a weak r 6 attracting part (van der 
Waals force); the spheres have two different radii to avoid 
crystallization. The mixture is a so called fragile glass 
former that has a glass transition at around T = 0.46 in 
the appropriate LJ-units. 



Fig. 2. Generalized scattering function C(t w , tw +t) as a function 
of the time for the quenches from 7 000 K (top) and 8 000 K 
(bottom) to 3 000 K and various selected waiting times t w . 
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In recent simulations of this model, Kob and Barrattf 
found a power law for the potential energy with an expo- 
nent of —0.144, which is of the same order of magnitude 
as our result for silica. We note, however, that in sim- 
ulations for a soft sphere modelJIP thus with a potential 
similar to the Lennard- Jones one, but using a Monte 
Carlo method instead of molecular dynamics, the differ- 
ent value 0.7 has been found for the exponent; Kob and 
Barrat mention the possibility that the disagreement is 
not caused mainly by the comparatively small differences 
of the models, but by an influence of the microscopic dy- 
namics on the aging process. The fact that we have ob- 
tained a result similar to theirs by using the same kind of 
simulation, but a different potential, might be a further 
indication for the validity of this assumption. 



3. 2 Scattering junction 

Figure 2 shows two sample plots for the time develop- 
ment of the scattering function. A longer waiting time 
i w means that the system has already had more time to 
relax after the quench than in the case of a short one; the 
speed of the decorrelation therefore decreases with grow- 
ing t w , so that the curves for larger values of t w decay 
more slowly than those for short waiting times. 

For short times the particles move ballistically, for long 
times they show a diffusive motion due to collisions with 
each other; particularly for low temperatures a plateau 
becomes visible between these two ranges, which is symp- 
tomatic for the situation that particles are trapped in 
"cages" formed by their neighbors and need a compara- 
tively long time to escape from these cages. The temper- 
atures in our simulations were not low enough to make 
this behavior visible very clearly, but the tendency to 
form a plateau is recognizable in the uppermost curves 
of figure 2. 



Fig. 3. Data for the same quenches as in figure 2, rescaled to have 
all curves fall together at C = 0.5. 
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In passing we mention that the curves for longer wait- 
ing times in figure 2 also show another interesting prop- 
erty, a dip before the plateau or the final decay is reached, 
at a time of around 0.2 ps. This dip is related to the so- 
called boson peak, a dynamical feature at around 1 THz 
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(10 12 sec -1 ) in the spectra of many strong glass formers; 
several theoretical approaches have been proposed to ex- 
plain this peak, but until now no consensus has been 
reached. 

A natural scaling ansatz for two-time correlation func- 
tions is 

C(t w , t„ + t) = / s tat(t) + 5agc(tAr) , 

where the first term is the stationary part 
/stat (t) = lim C (t w , i w + i) 

and <7age(^) is a scaling function for the aging part, which 
depends only on the ratio between the time t and an ef- 
fective relaxation time t r that varies monotonically with 
the waiting time t w . In simple coarsening systems the 
off-equilibrium properties are governed by a single length 
scale, the domain size or correlation length £ that grows 
with the waiting time as £ cx i^, and t r is simply pro- 
portional to t w . However, in_aging experiments on glassy 
systems like polymer glassescP and various spin glassesH 
it was found that t r oc describes the data reasonably 
well when the exponent a is used as a fit parameter (usu- 
ally a is close to but significantly different from 1). 

To examine whether a similar statement can be made 
for silica, we defined the relaxation time i r of the scat- 
tering function as the time necessary to reach the value 
C = 0.5 and checked the scaling ansatz by plotting C 
over t/t r . Figure 3 shows the results for two of the sim- 
ulated quenches: except for the case of small waiting 
times, the curves fall together well for t/t x > 1, confirm- 
ing the ansatz. 

Fig. 4. Relaxation time t T as a function of the waiting time i w for 
the quenches with initial temperature 6 000 K. The dotted lines 
represent the approximation by a power law for long waiting 
times (from 0.8 ps upwards) with the exponent a as given in 
table IV. 
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Plotting the relaxation time as a function of the wait- 
ing time reveals that, again with the exception of small 
waiting times, they can be approximated by a power law 
as well, as figure 4 demonstrates for the quenches with an 
initial temperature of 6 000 K. The exponents are listed 



in table IV; the average value turns out to be a = 0.89. 



Table IV. Estimated exponents a for the waiting time depen- 
dence of the relaxation time for the various combinations of ini- 
tial (rows) and final (columns) temperatures. (Due to an error 
when the simulations were run, not enough data could be gath- 
ered to give reliable results for the quench from 8 000 to 3 200 K.) 





2 700 K 


3 000 K 


3 200 K 


6 000 K 


0.88 


0.88 


0.85 


7000 K 


0.86 


0.95 


0.77 


8 000 K 


0.97 


0.96 





Again we compare our data to the one obtained, in 
simulations of the LJ model. Recent results thered'EF 
present a somewhat controversial picture of the behavior 
of the scattering function, which we summarize in figure 
5. 

Fig. 5. Scaling plot of the scattering function (for q = 7.2) for a 
Lennard-Jones system with 8000 particles quenched from T = 5 
to T = 0.35 (glass transition at T = 0.46). Top: The data 
are scaled according to a t/t r aging scenario, as in figure 3 for 
the silica system. The data collapse is not very good. Bottom: 
The data-pxe scaled according to a m((t+t w )/T)/ ln(tu,/T) aging 
scenario.^ the data collapse appears to be better. 
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[ ln((t+t w )/x) / ln(t w /T) ] - 1 



Here the data collapse is slightly better when using an 
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activated dynamics ansatzE 
+ t) =/ rt »t(t) 



/ln((* + <w)/ 



V ln(t w /r) 



although for higher temperatures the t/t r scaling sce- 
nario also works fine with t T ~ tf£ and a 0.88,tf sur- 
prisingly similar to our result for silica reported above. 

To conclude we reported results of molecular dynamics 
simulations of a model for a strong glass former, amor- 
phous silica, and showed 1) that the system is clearly 
aging on time scales of pico-seconds and 2) that the off- 
equilibrium or aging part of the generalized structure 
function obeys a t/t T scaling scenario with the effective 
relaxation time (or effective age) obeying t T oc with 
a « 0.88. These results are in good agreement with 
those reported recently for a fragile glass former, a bi- 
nary Lennard- Jones mixture. 
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